typedef float ft;
const int my_L = 128; // maximum 1024
const int my_M = 128;
const int my_N = 128;

// GPU kernel to add the elements of two arrays
template <typename T>
__global__ void gpu_version2(const T * __restrict__ input, T * __restrict__ output, const T * __restrict__ matrix, const int L, const int M, const int N){
  // parallelize threadIdx.x over vector length, and blockIdx.x across k (N)
  __shared__ T smem[my_L];
  int idx = threadIdx.x;
  int k = blockIdx.x;
    T v1 = 0;
    for (int i = 0; i < M; i++)
      v1 += input[k*M*L+idx*M+i];
    v1 /= M;
    for (int i = 0; i < L; i++){
      __syncthreads();
      smem[threadIdx.x] = v1 * matrix[i*L+idx];
      for (int s = blockDim.x>>1; s > 0; s>>=1){
        __syncthreads(); 
	if (threadIdx.x < s) smem[threadIdx.x] += smem[threadIdx.x+s];}
      if (!threadIdx.x) output[k+i*N] = smem[0];}
}

//template <typename T>
void test_kernel_launch(
    float *d_input, float *d_output, float *d_matrix,int L,int M, int N,
    cudaStream_t stream)
{
    float _min = 0.0f;              
    float _max = 1.0f;             

    int block_size = L;
    int grid_size = N;

  
    gpu_version2<<<grid_size, block_size, 0, stream>>>(
        d_input, d_output, d_matrix, L, M, N);
}
